
conf = omni_load_conf()
s = omni_read_cat(conf.cat)
ancil = omni_read_conffile('./conffiles/ancillary.conf')
local = omni_read_conffile('./conffiles/local_layout.conf')
v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart
restore,'./local/'+conf.survey+'_velocities.sav',/ver


vind = where(v.grs.vlsr GE 30. AND v.vlsr GE 30. AND v.grs.flag EQ 2 AND $
             v.l GE 25 AND v.l LE 29, nv)
print,nv

vind = where( v.cnum EQ 6289, nv)
vind = where( v.cnum EQ 5341, nv)
print,nv

vs = v[vind]
sfull = s

;;========================================
;; Loop to find an interesting source

FOR jj=0,nv-1 DO BEGIN
   
   ;; dl = abs(s.glon - 32.10)
   ;; db = abs(s.glat + 0.07)
   print,'VCNUM: ',vs[jj].cnum
   ind = where(sfull.cnum EQ vs[jj].cnum, ni)
   print,ni
   
   s = sfull[ind]
   
   IF n_elements(grs) EQ 0 THEN $
      restore,'./local/'+conf.survey+'_grs_spectra_r46.sav',/ver
   g = grs[ind]
   
   bgps = readfits('/d3/DIST_LOCAL/postage/'+conf.survey+'_data'+$
                   string(s.cnum,format="(I04)")+'.fits',hd)
   label = readfits('/d3/DIST_LOCAL/postage/'+conf.survey+'_label'+$
                    string(s.cnum,format="(I04)")+'.fits',lhd)
   
   ;; Find where GRS spectrum exists
   ig = where(total(grs.onspec,1) NE 0.,ng)
   print,ng
   
   ;; Find where HCO+ or NH3 velocities exist
   ih = where(v.mol[0].nspec NE 0, nh)
   in = where(v.mol[3].nspec NE 0, nn)
   
   print,nh,nn
   
   
   print,s.name
   
   
   ;;**********************************************************************
   ;; PLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOTPLOT
   myps,'./grs_paper/grs_rindex_ex_A.eps',xsize=3.1,ysize=3.0,thick=1.5
   
   oncol = 'blu4'
   offcol = 'deep pink'
   
   ;;=============================================
   ;; PLOT ME UP!
   print,!x.omargin,!y.omargin
   !x.omargin = [+0.5,-2.4]
   !y.omargin = [-0.75,-1.5]
   multiplot,[1,1],/doxaxis,/doyaxis
   
   csf = 0.55
   
   undefine,pr
   pr = [-0.10,1.40]
   ;;pr = [-0.13,0.60]
   plotmap,bgps,hd,ct=3,XC=xc, YC=yc,charsize=1.0*csf,range=pr,axcol='bisque',$
           xtickinterval=0.02,ytickinterval=0.02
   print,'RANGE: ',pr
   
   cgContour,label=0,label,xc,yc,/over,color=oncol,c_thick=15,level=[0.5]
   cgContour,label=0,g.border,xc,yc,/over,color=offcol,c_thick=5,level=[0.5]

   myps,/done,/mp
   
   
   
   myps,'./grs_paper/grs_rindex_ex_B.eps',xsize=3.1,ysize=3.0,thick=1.5
   multiplot,[1,1],/doxaxis,/doyaxis

   
   ;;===============================
   ;; The Spectra
   
   v_std = findgen(conf.nvbin)*conf.deltav + conf.vstart
   
   cgPlot,v_std,g.onspec,color=oncol,charsize=1.0*csf,xr=[-5,135],$[-5,115],
          /xst,thick=2,$
          xtit='V!dLSR!n  [km s!u-1!n]',ytit='!u13!nCO(1-0)  T!dA!n  [K]',$
          /nodata,yr=[0,2.8]
   cgOplot,v_std,g.bdrspec,color=offcol,thick=3
   cgOplot,v_std,g.onspec,color=oncol,thick=3
   cgOplot,v_std,g.spectrum,thick=6
   
   vline,/h,ancil.grs_ta,color='blk5'
   
   ;;   vline,v[ind].grs.vlsr,color='goldenrod'
   vline,v[ind].vlsr,color='green',thick=3
   print,v[ind].grs.vlsr,v[ind].vlsr
   
   v1ind = where( abs(v_std - 40.) LE 5.)
   v2ind = where( abs(v_std - 83.) LE 5.)
   print,max(g.spectrum[v1ind],jj1),max(g.spectrum[v2ind],jj2),$
         max(g.spectrum[v1ind])/max(g.spectrum[v2ind])
   print,v_std[v1ind[jj1]],v_std[v2ind[jj2]]
   
   al_legend,/top,/right,linsize=2.5,linestyle=[0,0,0,2],$
             color=[oncol,offcol,'black','green'],$
             ['Source Spectrum','Border Spectrum',$
              'On-Off','V!dLSR!n from HCO!u+!n(3-2)'],$
             charsize=0.9*csf,box=0,thick=[3,3,6,3]
   
   cgAxis,xaxis=0,/xst,xtickformat='blank_axis'
   cgAxis,xaxis=1,/xst,xtickformat='blank_axis'
   
   ;; multiplot,/doyaxis,/doxaxis
   
   ;; ;;===============================
   ;; ;; Velocity-Integrated GRS postage stamp
   
   ;; ellstr = string(round(v[ind].l),format="(I0)")
   ;; readcol,local.grs,grsfn,format='a'
   ;; grsi = where(strmatch(grsfn,'*'+ellstr+'*'),ngi)
   ;; IF ~exist(gcube) THEN gcube = readfits(grsfn[grsi],ghd)
   ;; rdhd,ghd,s=h
   ;; dv = median(h.v[1:*] - h.v[0:*])
   
   ;; veli = where(h.v GE v[ind].vlsr - 2*v[ind].lw AND $
   ;;              h.v LE v[ind].vlsr + 2*v[ind].lw, nvel)
   ;; print,v[ind].vlsr,v[ind].lw,nvel,dv
   
   ;; onmap = total(gcube[*,*,veli],3) * dv
   ;; help,gcube,onmap
   
   ;; mhd = ghd
   ;; sxaddpar,mhd,'NAXIS',2
   ;; sxdelpar,mhd,'NAXIS3'
   ;; sxaddpar,mhd,'EQUINOX',2000.
   ;; print,mhd
   
   ;; hastrom,onmap,mhd,pson,pshd,hd
   
   ;; onpr = set_plot_range(pson)
   ;; onpr[0] = 0.
   ;; plotmap,pson,pshd,ct=13,range=onpr,charsize=1.0*csf,ytickinterval=0.04,$
   ;;         xtickinterval=0.04,title='On'
      
   ;; cgContour,label=0,label,xc,yc,/over,color=oncol,c_thick=15,level=[0.5]
   ;; cgContour,label=0,g.border,xc,yc,/over,color=offcol,c_thick=5,level=[0.5]

   ;; multiplot,/doyaxis,/doxaxis
   
   ;; ;;===============================
   ;; ;; "On-Off" GRS postage stamp
   ;; veli = where( (h.v GE v[ind].vlsr - 4*v[ind].lw AND $
   ;;                h.v LE v[ind].vlsr - 2*v[ind].lw) OR $
   ;;               (h.v GE v[ind].vlsr + 2*v[ind].lw AND $
   ;;                h.v LE v[ind].vlsr + 4*v[ind].lw), nvel)
   ;; print,v[ind].vlsr,v[ind].lw,nvel,dv
   
   ;; offmap = total(gcube[*,*,veli],3) * dv
   ;; help,gcube,offmap
   
   ;; hastrom,offmap,mhd,psoff,pshd,hd
   
   ;; plotmap,pson-psoff,pshd,ct=13,range=onpr,charsize=1.0*csf,$
   ;;         ytickinterval=0.04,xtickinterval=0.04,title='On-Off'
   
   ;; cgContour,label=0,label,xc,yc,/over,color=oncol,c_thick=15,level=[0.5]
   ;; cgContour,label=0,g.border,xc,yc,/over,color=offcol,c_thick=5,level=[0.5]

   ;; print,onpr
   
   ;; cgColorbar,position=[0.3,0.03,0.8,0.05],charsize=0.6,range=onpr,$
   ;;            title='Integrated !u13!nCO Intensity  [K km s!u-1!n]'
   
   
   myps,/done,/mp

   ;; multiplot,/doyaxis
   
   ;; ;;===================================================
   ;; ;; HCOP Comparison panel
   
   
   ;; ;; Get list of GRS-HCOP compatabilities
   ;; iginh = where_array(ig,ih,ngh)
   ;; ihing = where_array(ih,ig,nhg)
   ;; print,ngh,nhg
   ;; ihp = ih[iginh]
   ;; igp = ig[ihing]
   ;; ;; Update ihp to remove grs.flag = 0
   ;; igz = where(v[ihp].grs.flag NE 0 AND v[ihp].grs.flag NE 3, ngz)
   ;; print,ngz
   ;; ihp = ihp[igz]
   
   ;; cgPlot,xr=[-50,150],yr=[-50,150],v[ihp].vlsr,v[ihp].grs.vlsr,psym=16,$
   ;;        symsize=0.4,/xst,$
   ;;        ytit='On-Off !u13!nCO V!dLSR!n [km s!u−1!n]',charsize=1.0*csf,$
   ;;        xtit='HCO!u+!n  V!dLSR!n  [km s!u-1!n]'
   
   ;; vline,-5,color='blk5'
   
   ;; ineg = where(v[ihp].mol[0].vlsr LT -5, nneg)
   ;; cgOplot,v[ihp[ineg]].vlsr,v[ihp[ineg]].grs.vlsr,$
   ;;         psym=16,symsize=0.41,color='Deep Pink'
   
   ;; ;; Compute matching to w/in 5 km/s
   ;; ikms = where(abs(v[ihp].vlsr-v[ihp].grs.vlsr) LE 5., nkms)
   ;; al_legend,/top,/right,box=0,charsize=0.9*csf,$
   ;;           [string(float(nkms)/n_elements(ihp)*100.,format="(F0.1)")+$
   ;;            '% with '+cgSymbol('Delta')+'V!dLSR!n '+cgSymbol('leq')+$
   ;;            ' 5 km s!u-1!n']
   ;; al_legend,/bottom,/right,box=0,['N = '+string(ngz,format="(I0)")]
   
   ;; ;; al_legend,/top,/left,box=0,['(c)'],charsize=1.2

   ;; myps,/done,/mp


   ;; res = get_kbrd()
ENDFOR

END
